

######################################################

  Table 3. Summary of Example 1 Test Statistics

#####################################################

library(foreign)
library(lavaan)
library(boot)

mw03<-read.dta("/Users/bv533/Dropbox/LM-Test/new empirical example/stony_brook/mw03_survey2.dta")


dim(mw03)


names(mw03)


model<-"ni=~q22 + q20 + q23 + q24 
        sp=~q28 + q26 
        up=~q22 + q40 + q27 + q39 + q37 + q38 
"

newdata2 <- lapply(mw03,as.numeric)
newdata2<-data.frame(newdata2)
newdata2

f01<-cfa(model, data=newdata2, ordered=c("q22", "q20", "q23", "q24", "q28", "q26", 
                                         "q22", "q40", "q27", "q38", "q39", "q37"))
f01



modindices(f01, sort = TRUE)[1:10,]

newpar=' 

sp=~q27
ni=~q27
q27~~q38
q37~~q38
ni=~q40
q40~~q38
q40~~q39
sp=~q37
q37~~q39
q27~~q39'


## This prints out Univariate LM test and p-value in Table 3. 

lavTestScore(f01, add=newpar)




############ Bootstrap LM Test

newpar=' 
sp=~q27
ni=~q27
q27~~q38
q37~~q38
ni=~q40
q40~~q38
q40~~q39
sp=~q37
q37~~q39
q27~~q39'



set.seed(112)

Data<-newdata2
bin.1000 <- matrix(NA, nrow = 500, ncol = 2)
for (i in 1:500) {
  
  boot.idx <- sample.int(nrow(Data), replace = TRUE)
  Data.boot <- unique(Data[boot.idx,])
  
  f02 <- cfa(model = model, data = Data.boot, ordered=c("q22", "q20", "q23", "q24", "q28", "q26", 
                                                        "q22", "q40", "q27", "q38", "q39", "q37"))
  
  bin.1000[i, 1] <- lavTestScore(f02, add = newpar)$uni[1, 4]
  bin.1000[i, 2] <- lavTestScore(f02, add = newpar)$uni[1, 6]
}

t1 <- mean(bin.1000[, 1], na.rm = TRUE)
t2 <- mean(bin.1000[, 2], na.rm = TRUE)

t1
t2



####. bootstrap LM chis-square test statistic and p-value will print here



############ Bootstrap Wald Test


w_model<-"ni=~q22 + q20 + q23 + q24 + b2*q27 + b5*q40
        sp=~q28 + q26 + b1*q27
        up=~q22 + q40 + q27 + q39 + q37 + q38 

q27~~b3*q38
q37~~b4*q38



"

con='

b1==0

'

set.seed(115)

Data<-newdata2
bin.1000 <- matrix(NA, nrow = 500, ncol = 2)
for (i in 1:500) {
  
  boot.idx <- sample.int(nrow(Data), replace = TRUE)
  Data.boot<-unique(Data[boot.idx,]) 
 

  fit2 <- cfa(model = w_model, data = Data.boot, ordered=c("q22", "q20", "q23", "q24", "q28", "q26", 
                                                         "q40", "q27", "q38", "q39", "q37")) 
  
  bin.1000[i, 1] <- lavTestWald(fit2, constraints = con)$stat[1]
  bin.1000[i, 2] <- lavTestWald(fit2, constraints = con)$p.value[1]
  
}

w1 <- mean(bin.1000[, 1], na.rm = TRUE)
w2 <- mean(bin.1000[, 2], na.rm = TRUE)

w1
w2







############################.  This is for Table 4 

new_model<-"ni=~q22 + q20 + q23 + q24 + q27
        sp=~q28 + q26
        up=~q22 + q40 + q27 + q39 + q37 + q38


"


f03 <- cfa(model = new_model, data = newdata2, ordered=c("q22", "q20", "q23", "q24", "q28", "q26", 
                                                          "q40", "q27", "q38", "q39", "q37"))
f03
fitMeasures(f03, c("chisq", "df", "pvalue", "nfi", "tli", "cfi", "rmsea"))



##########. The End

##########. All results match the manuscript, 12/6/24
